############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)
base::library(haven)
base::library(nnet)
base::library(ggeffects)

load(here("Data", "OtherQuestions.RData"))

summary(m1 <- lm(law_enforce ~ death_rate*party_id3 + unemp_rate + defl_pcincome +
                    log_pop + rural_urban + defl_pcincome +
                    gender + black + hispanic + asian + other_race + 
                    log_age + college + income + year_elec, 
                  weights = weight_cumulative, data = ces_data))

summary(m2 <- lm(healthcare ~ death_rate*party_id3 + unemp_rate + defl_pcincome +
                   log_pop + rural_urban + defl_pcincome +
                   gender + black + hispanic + asian + other_race + 
                   log_age + college + income + year_elec, 
                 weights = weight_cumulative, data = ces_data))

summary(m3 <- glm(more_borderpatrols ~ death_rate*party_id3 + unemp_rate + defl_pcincome +
                    log_pop + rural_urban + defl_pcincome +
                    gender + black + hispanic + asian + other_race + 
                    log_age + college + income + year_elec, 
                  weights = weight_cumulative, 
                  family = binomial(link = "logit"), data = ces_data))


model_border <- m3 %>% 
  ggeffect(terms = c("death_rate[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120]", 
                     "party_id3"), 
           ci.lvl = 0.90)
glimpse(model_border)
table(model_border$response.level)

model_border <- data.frame(model_border)
glimpse(model_border)

#What do you think the U.S. government should do about immigration? Increase the number
#of border patrols on the US-Mexican border.
fig_border <- model_border %>% 
  #filter(response.level == "Rep") %>% 
  mutate(group = factor(group,
                        levels = c("Democrat", 
                                   "Independent",
                                   "Republican"),
                        labels = c("Democrats", 
                                   "Independents",
                                   "Republicans"))) %>% 
  ggplot(aes(shape = group, color = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85#, position = position_dodge(width = .5),
                  #color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability",
       title = NULL) +
  ggtitle("C. Increase the Number of Border Patrols") +
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c("blue", "gray40", "red")) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  scale_y_continuous(breaks = c(0.2, 0.4, 0.6, 0.8, 1), lim = c(0.2, 1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))
fig_border

model_lawenf <- m1 %>% 
  ggeffect(terms = c("death_rate[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120]", 
                     "party_id3"), 
           ci.lvl = 0.90)
class(model_lawenf)
table(model_lawenf$response.level)

model_lawenf <- data.frame(model_lawenf)
glimpse(model_lawenf)
class(model_lawenf)

fig_law <- model_lawenf %>% 
  #filter(response.level == "Rep") %>% 
  mutate(group = factor(group,
                        levels = c("Democrat", 
                                   "Independent",
                                   "Republican"),
                        labels = c("Democrats", 
                                   "Independents",
                                   "Republicans"))) %>% 
  ggplot(aes(shape = group, color = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85#, position = position_dodge(width = .5),
                  #color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Support",
       title = NULL) +
  ggtitle("A. Increase Spending on Law Enforcement") +
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Support:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c("blue", "gray40", "red")) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  scale_y_continuous(breaks = c(1.8, 2.2, 2.6, 3, 3.4),
                     lim = c(1.7, 3.5)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))
fig_law


model_healthcare <- m2 %>% 
  ggeffect(terms = c("death_rate[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120]", 
                     "party_id3"), 
           ci.lvl = 0.90)
class(model_healthcare)

model_healthcare <- data.frame(model_healthcare)
glimpse(model_healthcare)
class(model_healthcare)

fig_hc <- model_healthcare %>% 
  #filter(response.level == "Rep") %>% 
  mutate(group = factor(group,
                        levels = c("Democrat", 
                                   "Independent",
                                   "Republican"),
                        labels = c("Democrats", 
                                   "Independents",
                                   "Republicans"))) %>% 
  ggplot(aes(shape = group, color = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85#, position = position_dodge(width = .5),
                  #color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Support",
       title = NULL) +
  ggtitle("B. Increase Spending on Healthcare") +
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c("blue", "gray40", "red")) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  scale_y_continuous(breaks = c(1.8, 2.2, 2.6, 3, 3.4), 
                     lim = c(1.7, 3.5)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))
fig_hc

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(fig_law, fig_border, fig_hc, cols = 2)


summary(m2)
stargazer(m1, m2, m3,
          type = "latex",
          title = "Statistical Models --- Figure 6",
          dep.var.labels = c("Law Enforcement", 
                             "Healthcare", 
                             "Border Patrols"), 
          column.labels = c("Panel A", 
                            "Panel B", 
                            "Panel C"),
          no.space = FALSE, single.row = T,
          covariate.labels = c("Overdose Death Rate",
                               "PID3 - Independent",
                               "PID3 - Republican",
                               "Local Unemployment Rate",
                               "Per Capital Income (County",
                               "Total Population (Log)",
                               "Gender - Woman",
                               "Race - Black",
                               "Race - Hispanic",
                               "Race - Asian",
                               "Race - Other",
                               "Age (Log)",
                               "College Educated",
                               "Household Income",
                               "Death Rate X PID3 - Independent",
                               "Death Rate X PID3 - Republican",
                               "Constant"),
          omit = c("rural_urban", "year_elec"),
          omit.labels = c("Rural-Urban Code",
                          "State Fixed Effects"))

